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ABSTRACT 

Magnetohydrodynamics (MHD) is invoked to address turbulent fluctuations in a variety of astro- 
physical systems. MHD turbulence in nature is often anisotropic and imbalanced, in that Alfvenic 
fluctuations moving in opposite directions along the background magnetic field carry unequal energies. 
This work formulates specific requirements for effective numerical simulations of strong imbalanced 
MHD turbulence with a guide field Bp. High-resolution simulations are then performed and they 
suggest that the spectra of the counter-propagating Alfven modes do not differ from the balanced 
case, while their amplitudes and the corresponding rates of energy cascades are significantly affected 
by the imbalance. It is further proposed that the stronger the imbalance the larger the magnetic 
Reynolds number that is required in numerical simulations in order to correctly reproduce the turbu- 
lence spectrum. This may explain current discrepancies among numerical simulations and observations 
of imbalanced MHD turbulence. 
Subject headings: 



1. INTRODUCTION 

Magnetohydrodynamic (MHD) turbulence plays a 
critical role in theoretical modeling of the spectra 
of velocity and magnetic fluctuat ions in the Solar 
Wind (e.g.. IColemanlfl96l. [1 968: Be lcher fc DavislHOTl 
IMarschl [19911 : [Goldstein etal.iil995i ) as well as electron 
density fluctuations in the interstellar medium (e.g., 
[Armstrong eFallfTggsl : fLithwick fc Go ldreich"2001'). Fol- 
lowing the pionee ring works of llrcSmikov ()1963i) and 
iKraichnanl ()1965[ ). and t he introduction of t h e con - 
cept of critical balance bv iGoldreich fc SridhaH (|1995f ). 
progress in MHD turbulence was mostly due to phe- 
nomenological models, closure theories, and weak tur- 
bulence asymptotic models. In recent years, constantly 
increasing computing power has put within our reach 
the ability to simulate universal inertial range spectra 
of strong MHD tur bulenc e (e.g.. iCho fc Vishniac 2000; 
Maron fc Goldreich' '2001: Haugcn et al. 2003; Biskamo 
2003; Miiller fc Graooin 2005; Mason ct al. 2006, 200&; 
Mininni fc Pou auct 2007; Beresnvak fc LazarianI 120081 : 
Perez fc Boldvrev 2008. 200^ These simulations have 



motivated new phenomenological models and have re- 
vived the vigorous interest in the fundamentals of MHD 
turbulence. 

As an important result, simulations have shown that 
strong MHD turbulence is locally imbalanced, that is, it 
spontaneously develops correlated regions of positive and 
negative cross-helicity where the energy in Alfven waves 
propagating along and against the magnetic field are not 
equal, irrespective of the total amount of cross-heli city in 
the sy stem, see Perez fc Boldyrev (2009); Boldvre v et al.l 
(|2009f l. In each of these regions both energy and cross 
helicity are subject to a nonlinear cascade from large to 
small scales, consequently having an effect in the overall 
energy spectrum. 

Such an imbalance is clearly present in the solar wind, 
where velocity and magnetic fluctuations show high cor- 



relations of a preferred sign, that is, the normalized cross 
helicity ac = {v ■ h) / E = Hc/E is predominantly close 
to unity. Here E is the average total energy (kinetic plus 
magnetic), He is the cross helicity and the brackets de- 
note a suitable ensemble average. The preferred positive 
sign of (Tc indicates that there is more energy in Alfven 
waves propagating outwards from the Sun than propa- 
gating inwards. 

Several phenomenological models have been pro- 
posed to addres s stro n g ini balanced MHD tur- 
bulence ( Lithwick el_al.l 1200^ 
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Beresnvak' fc Lazarian (2008); P erez fc Boldvrev! ([2009); 
Podesta fc Bhattachariee (,2009i) ). some with support 
from numerical simulations and observations. However, 
these works have led to conflicting predictions. For in- 
stance, the theory by Lithwick et al. concludes that in 
imbalanced regions the Elsasser spectra have the same 

scalings E^(kx) cx E~{k±) oc kj_ . The theory by 
Chandran proposes that the spectra of E^{k±) and 
E~{k±) are different depending of the degree of imbal- 
ance. The theory by Beresnyak fc Lazarian also sug- 
gests the different spectra for E^{k±) and E~{k±). Fi- 
nally, the analysis by Perez fc Boldyrev and Podesta 
fc Bhattacharjee finds that the spectra of E~^{k±) and 
E^{k±) have different amplitudes but the same scalings 

E+{kj_) (xE-{k±) cx k^^^"^. 

On the numerical side, anisotropy and imbalance sig- 
nificantly complicate simulations of MHD turbulence. In 
the presence of a strong guide field, the turbulent fluc- 
tuations become elongated in the direction of the field. 
This increases the variety of numerical settings that can 
be used to drive such turbulence and complicates com- 
parison of numerical simulations carried out by differ- 
ent groups. For example, numerical settings can include 
weak or strong guide fields; different levels of imbalance; 
cubic or rectangular simulation boxes; short or long-time 
correlated large-scale forces; frozen large-scale modes; 
forcing that excites only v, b or both; volume forcing 
or forcing at the boundaries of the domain; different res- 
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olutions along field-perpendicular and field-parallel di- 
rections, etc. 

To understand how different numerical simulations can 
be compared to each other, one needs to understand the 
conditions that should be met to ensure that simulations 
reproduce the range of scales where turbulence is univer- 
sal. For instance, one has to be sure that the forcing and 
the dissipation do not spoil the inertia l interval. The first 
steps in this direction w ere made in (jMason et al.ll2006l : 
iPerez &: BoldvrevI l2008f ) for the balanced case. It was 
found that the inertial range scaling can be sensitive to 
the way in which the turbulence is forced, not because of 
a lack of universality, but because the inertial range has 
a limited extent. 

In the present paper we formulate the conditions that 
need to be satisfied in order to correctly simulate strong 
imbalanced MHD turbulence. We then propose a univer- 
sal numerical setting, based on the reduced MHD model, 
which allows one to find the spectra of both balanced and 
imbalanced MHD turbulence. We found that imbalance 
significantly reduces the inertial interval in numerical 
simulations. The numerical resolution required to pro- 
duce a large inertial interval strongly increases with the 
amount of imbalance. An imbalance essentially stronger 
than 7^ = ^ 10 cannot be accessed with 

currently available supercomputers, which significantly 
limits the applicability of present-day numerical simula- 
tions to practical cases. 

2. MODEL EQUATIONS 

In the presence of a guide field Bp (say in the z di- 
rection), the MHD equations describing the evolution of 
magnetic and velocity fluctuations, b(x, t) and v(x, t) 
can be written in terms of the so-called Elsasser vari- 
ables, = V ± b: 

9tz± T (Va • V)z± + (zT • V)z± = -VP + f , (1) 

where = Bo/\/47rp is the Alfven velocity, p is the 
fluid density, P is the total pressure that is determined 
from the incompressibility condition, V • z^ = 0, f is 
large-scale forcing and we neglect small viscosity and re- 
sistivity. In the limit of weak fluctuations, the incom- 
pressible MHD system describes non-interacting Alfven 
waves with dispersion relation uj^(k.) — ±fc||VA. These 
waves can have two types of polarizations, the shear- 
Alfven one (z^) and the pseudo- Alfven one (zp) given 
in Fourier space by eg = x k/fc_L and ep = k x eg/fc, 
where e|| is a unit vector in z direction. Strong MHD 
turbulence is dominated by fluctuations with k± ^ 
iGoldreich fc Sridhail (|1995f ) argued that since for large 
k± the polarization of the pseudo- Alfven fluctuations is 
almost parallel to the guide field, such fluctuations are 
coupled only to field-parallel gradients, which are small 
since fc|| <C k±. Therefore, the pseudo- Alfven modes do 
not play a dynamically essential role in the turbulent 
cascade. We can remove the pseudo-Alfven modes by 
setting Z|j^ = in equations ([T|) to obtain 

atZ±T(VA-V)z±-f (iT.V)i± = -ViP+iv2z±, (2) 

K 

where R denotes the Reynolds numbers (discussed be- 
low). In this system, the fluctuating fields have only 
two vector components, = {z^,z^,0}, but depend 



on all three spatial coordinates. System ^ is equivalent 
to the Reduced MH D model, originally develope d for 
tokamak plasm as by iKadomtsev fc Pogutsd (|1974f ) and 
■Straus^ HFM ). 

3. NUMERICAL STRATEGY 

Model ([2]) describing the nonlinear interactions of 
Shear- Alfven modes is a valuable tool in theoretical and 
numerical studies of incompressible MHD turbulence. 
Based on this model, we now discuss the conditions that 
should be satisfied in order to correctly simulate strong 
MHD turbulence, both balanced and imbalanced. 

3.1. Computation Domain 

In the presence of a strong gu ide field the turbu- 
lence is anisotropic. According to IGoldreich fc Sridhaii 
(1995), deep in the inertial range turbulence should be- 
come strong and the fluctuations should satisfy the crit- 
ical balance condition fcyPo ~ k±bx, where fcy ~ l/l and 
k± 1/A are the fleld-parallel and field-perpendicular 
wave vectors associated to an anisotropic eddy of paral- 
lel size I and perpendicular size A, respectively. Let us 
define the nonlinear interaction strength parameter 

X^{k^bx)/ik\\Bo); (3) 

the critical balance condition then implies % ~ 1. 

Simulations of incompressible turbulence are generally 
based on the Fourier pseudo-spectral method. In order 
to allow for the inertial interval to develop, turbulence 
is driven at the lowest resolvable wave numbers, and the 
energy dissipates at large wave numbers determined by 
the Reynolds numbers. For simulations on a cubic pe- 
riodic box of size L, the smallest wave-numbers along 
the field-parallel and field-perpendicular directions coin- 
cide, i.e., k± = /cj| = 27r/L. Therefore, driving at the 
low /c|| , k± results in an isotropic forcing and the non- 
linear strength parameter at the forcing scale becomes 
Xo = k±b\/k^\Bo ~ bx/Bo <C 1, which means that 
at least at the large scales, nonlinear interactions are 
weak. This would not be harmful if we had the resources 
to achieve arbitrarily high resolution, as the turbulence 
would proceed weakly until % ~ 1, and then would be- 
come strong. However, simulations generally produce a 
rather limited inertial range, so that the parameter x can 
hardly reach unity in s uch a setup. 

As pointed out by iMaron fc GoldreichI (1200 iP). and 
applied in recen t siin ulations iMason et al.l ^006) ; 
IPerez fc Boldvre^ ^00§), the best way to avoid this is 
to use an anisotropic domain such that at the forcing 
scale the parameter x is already of order unity, that is, 
the excited large-scale modes are already anisotropic and 
satisfy the critical balance condition. To achieve this we 
choose an elongated box x L||, so that the lowest 
field-perpendicular and field-parallel wave-numbers are 
k± — 2tt/L± and fc|| = 27r/Pj|, respectively. In this case, 
forcing at the lowest fc_L,fc|| leads to xo = L\\bxl L^^B^, 
which is of order unity provided that 

^ bx/Bo. (4) 

In this way, the turbulence is excited in a strong regime 
and the cascade proceeds down to smaller scales preserv- 
ing the critical balance condition. 
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3.2. Numerical Resolution 

At first sight, it appears that elongating the box along 
the z direction to match the elongation of the eddies 
should not change the number of grid points required 
in this direction compared to the number of points in 
the X and y directions. Fortunately, the number of 
points in the z direction can be reduced. This fol- 
lows from the fact that the turbulent spectrum de- 
clines quite slowly, as a power-law, in the direc- 
tion, while it drops sharply in the A;|j direction for 
A; 1 1 > A;", where a is a some positive power not exceed- 
ing 1 (ICho fc yishniacl[2000l: IM aron k Goldreichl[200l 
lOuehton et all I2004t IPerez k BoldvrevI 120081 12009^ . 
This qualitatively different spectral behavior in A;|| and 
fc_L directions allows one to reduce the numerical resolu- 
tion by a factor of 2 to 4 in the parallel direction, see 
Table [T] We checked that the restoration of the full res- 
olution in the z direction does not change the results, 
while significantly increases the computing costs. 

3.3. Periodic Boundary Conditions 

The spectral method assumes periodic boundary con- 
ditions in all spatial directions. The periodic bound- 
ary conditions in the z-direction may cause questions of 
whether the Alfven modes counter-propagating along a 
given magnetic field line will always interact only among 
themselves. The answer is no, since periodicity of the 
fluctuations does not imply periodicity of the magnetic 
field lines. Magnetic field lines are not periodic, and each 
given eddy interacts with many independent counter- 
propagating eddies. One can check numerically, that 
reducing the parallel box size L\\ below (U) somewhat 
spoils the spectrum at low wave numbers, e.g., as seen 
in forced simulations of ([Miiller & Grappin 2005), while 
increasing it beyond ^ does not change the results, but 
increases the computational cost [Mason & Cattaneo, un- 
published] . 

3.4. Reynolds numbers 

Probably the most significant limitation is imposed 
by the consideration of imbalanced turbulence. Indeed, 
in the imbalanced case, 7 — z'^ / z^ > 1, the formal 
Reynolds numbers corresponding to z+ and tT fields are 
essentially different. Therefore, the resolution require- 
ments increase with the amount of imbalance in order to 
produce large inertial ranges. Assume that the number of 
grid points in a field-perpendicular direction scales with 
the Reynolds number as ~ Re^ , where /3 = 2/3 or 3/4 
depending on the spectral slope (3/2 or 5/3). Then in- 
creasing the imbalance 7 by, say, a factor of 3 will require 
increasing the resolution by approximately a factor of 2. 
Noting that resolution of at least 1024 is required to sim- 
ulate the imbalance 7 ~ 2 (see below), we conclude that 
significantly stronger imbalance is not achievable with 
present day computing power. 

3.5. Random Forcing 

In this section we discuss the important aspects to be 
considered when choosing a particular forcing. We as- 
sume a random force f that has no component along z, 
it is solenoidal in the x — y plane and its Fourier coeffi- 
cients outside the range 1 < fc_L < 2, (27r/L||) < fc|| < 



Run 


Resolution 


Re 


u = 'q 








Al 


512^ X 256 


2400 


4.2 X 10" 


-i 


5 





A2 


10242 X 256 


6000 


1.7 X 10" 


-4 


5 





Bl 


256=* 


900 


1.1 X 10- 


-3 


10 


0.6 


B2 


512^ X 256 


2200 


4.6 X 10" 


-4 


10 


0.6 


B3 


10242 X 256 


5600 


1.8 X 10" 


-4 


10 


0.6 



TABLE 1 

Summary of simulations of strong balanced turbulence 
(al, a2) and strong imbalanced turbulence (bl, b2, b3). 

(27r/L||)n2 are zero, where Uz determines the width of 
the force spectrum in fc|| , and Lj_ = 2'k. The Fourier co- 
efficients inside that range are Gaussian random numbers 
with amplitude chosen so that the resulting rms velocity 
fluctuations are of order unity. The individual random 
values are refreshed independently at time intervals r. 
The parameter controls the degree to which the crit- 
ical balance condition is satisfied at the forcing scale. 
Note that we do not drive the fc|| = mode but allow it 
to be generated by nonlinear interactions. 

In contrast with incompressible hydrodynamic sys- 
tem, an incompressible MHD system can support Alfven 
waves. When the system is driven by a time depen- 
dent forcing, the most effectively driven modes are those 
resonating with the frequencies present in the forcing. 
Therefore, the spatial spectra of the large-scale velocity 
and magnetic fluctuations are not generally the same as 
the spatial spectrum of the force. Rather, they essen- 
tially depend on the both spatial and temporal spectra 
of the random forcing. Ultimately, it is the spectrum of 
the large-scale fluctuations, not the driving force, that 
should be controlled in numerical simulations. 

In the following we perform simulations with a short- 
time-correlated random forcing that drives the turbu- 
lence close to the critical balance condition at the large 
scales. The short correlation time t ensures a broadband 
frequency spectrum for the forcing, which allows high 
frequency Alfven waves to be excited, so that one can 
control the width of the field-parallel spatial spectrum of 
the fluctuations by choosing the number of driven field- 
parallel modes, Uz- A sort-time correlated Gaussian ran- 
dom force has another important advantage. It allows 
one to control the rates at which the energies of z+ and 
z" modes are injected. 

4. NUMERICAL RESULTS 

Table [1] summarizes five representative simulations 
that incorporate all the aspects discussed in section [3l 
The simulations produce consistent and physically mean- 
ingful results for a range of Reynolds numbers and de- 
grees of imbalance. Runs A are carried out for balanced 
turbulence, that is, Cc = 0. In runs B, cross helicity 
is injected at the forcing scale in such a way that CTc 
reaches a steady state of ^ 0.6. We use a short time- 
correlated forcing, compared to the Alfven time of the 
excited modes, so that the energy injection rates for both 
z"*" and z~ only depend on the variance of the imposed 
forcing, which is controlled in our simulations. In the 
imbalanced case, field-parallel box size is optimized to 
reach the critical balance at the large scales. Except for 
the Reynolds numbers, simulations Bl, B2, and B3 have 
the exact same parameters including the energy injection 
rates, e"*" and e~ . 

Since the background magnetic field must be strong. 
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we choose Bq = 5 in the Vrms units, of the discussion 
m ason et aLll2006f ). Time is normahzed to the large 
scale eddy turnover time tq = L±/2'iTVrms, where L± is 
the field-perpendicular box size. The Reynolds number 
is defined as Re = v^mg[L^/2'!T)/v and we have chosen 
the same value for the magnetic Reynolds number, Rm = 
brms{L±/2Tr)/r], denoting both by i? in ((21). In each run, 
the average is performed over about 100 large-scale-eddy 
turnover times. The results are presented in Fig. ((!]). 
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Fig. 1. — Spectra of the Elsasser fields in numerical simulations 
of MHD turbulence. Top two frames: balanced turbulence (runs 
Al, A2); bottom frame: imbalanced turbulence (runs B1-B3). 



5. DISCUSSION 

The advantage of our optimized setup (reduced MHD, 
elongated box, reduced resolution in z-direction) can 
be seen already in simulations of balanced turbulence, 
top frames in Fig. ([T]), runs Al and A2. The energy 

spectra approach E^{k±) oc fc, in good agreement 
with e arlier numerica l findin gs (e.g.. | M ar on fc Goldreic 



2m[ IHaugen et^al.1 IIBM IMiiller fc Grappinl 1200 
Mason et al.ll2006f) . while our simulations require con 



siderably less computational cost and produce slightly 
larger inertial intervals. Our most significant results 
are obtained for the imbalanced case. The bottom 
frame of Fig. [1] shows the spectra for three differ- 
ent Reynolds numbers (runs Bl to B3). It is ob- 
served that the spectra are pinned at the dissipa- 
tion scales, which supports the phenomenological pre- 



dictions by IGrappin etall (I1983D; IGaltier et all (| 2000[ ); 
ILithwick fc GoldreichI (120031 ) HChandraTT plOS). We also 
find that the large-scale parts of the both spectra are 
practically insensitive to the Reynolds numbers. These 
two important properties imply that as the Re numbers 
are further increased, the E spectra must become pro- 
gressively more parallel in the inertial interval. This is 
indeed seen in our numerical simulations (B1-B3). More- 
over, our numerical simulations indicate that both spec- 
tra approach the universal scaling of strong MHD turbu- 
lence E^{k±) cx k^^^^, while they have essentially dif- 
ferent amplitudes and correspond to essentially different 
energy ffuxes. 

The proposed numerical setup and the results of our 
simulations can serve as an important tool for resolving 
currently existing controversies regarding the spectra 
of imbalanced MHD turbulence. In particular, our 
results are consistent with the theories and observations 
predicting same sc ali ng — 3/2 for both Elsasser fields 
f Perez fc BoldvrevI I2009t IPodesta fc Bhattachariei 
2009). They are al s o bro adly consistent with the theory 
by Li thwick et al.l (|2007f ). predicting the same scaling 
—5/3 for both fields. The major difference between 
the two phenomenologies is in the phenomenon of 
dynamic alignment that is taken into account in the 
former models and is not considered in the latter one. 
Our numerical findings are less consistent with the 
models predicti ng different sca l ings for the Elsasser 
fields i?^, (e.g.. iChandranI [20081 : iBeresnvak fc LazarianI 
|2008() . It would be interesting to understand what 
assumptions of these models disagree with the numerics. 
Finally, our analysis m ay e xplain somewhat puzzling 
numerical findings by IBere snvak fc Lazarian (2009), 
who report different spectra for the E^ Elsasser fields, 
and the intersection of the spectra rather than pinning 
at the dissipation scale. According to our results, 
the explanation might lie in the fact that in these 
simulations the imbalance was extremely high, up to 
7^ — (z+)^/(2;~)^ ~ 1000, and, therefore, the universal 
regime of imbalanced MHD turbulence was not reached, 
see our analysis in Sec. 13.41 
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